Trophic ecology of Octopus vulgaris paralarvae along the Iberian Canary current eastern boundary upwelling system

Our knowledge of the diet of wild octopus paralarvae, Octopus vulgaris, is restricted to the first 2 weeks of its planktonic phase when they are selective hunters found near the coastline. These small paralarvae, bearing only three suckers per arm, are transported by oceanic currents from the coast towards offshore waters, where they complete the planktonic phase over 2 months. Here, we have investigated the trophic ecology of O. vulgaris paralarvae in two contrasting upwelling sub-regions of the Iberian Canary current (ICC) eastern boundary upwelling system and have evaluated dietary change as paralarvae develop (inferred by counting the number of suckers per arm, ranging from three to 15) along the coastal-oceanic gradient during their planktonic phase. Using high-throughput amplicon sequencing, we have characterised the diet of 100 paralarvae collected along the Northwest Iberian Peninsula (n = 65, three to five suckers per arm) and off the west coast of Morocco (n = 35, three to 15 suckers per arm), identifying up to 87 different prey species. The diet of paralarvae varied along the ICC, with crabs (53.4%), siphonophores (12.2%), copepods (12.3%), cnidarians (8.4%) and pteropods (3.7%) accounting for 90% of the variability detected off NW Iberian Peninsula, whereas off W Morocco, crabs (46.2%), copepods (23.1%), cnidarians (12.9%), krill (9.3%) and fishes (4.2%) explained 95.6% of the variability observed using frequency of observance (FOO%) data. Ontogenetic changes in the diet based on groups of paralarvae with similar numbers per arm were evidenced by the decreasing contribution of coastal meroplankton and an increase in oceanic holoplankton, including siphonophores, copepods, pteropods and krill. Trophic niche breadth values ranged from 0.06 to 0.67, with averaged values ranging from 0.23 to 0.33 (generalist = 1 and specialist = 0), suggesting that O. vulgaris paralarvae are selective predators through their ontogenetic transition between coastal and oceanic environments.


Results
Prey composition and abundance. Dissections of the digestive tract of 100 O. vulgaris paralarvae revealed that 91 had empty crops and nine stomachs with partially digested prey, none of which could be identified visually. Five paralarvae had no prey DNA detected; four were collected in Morocco (with three, four (n = 2) and six suckers per arm), and one was collected off the coast of Portugal with five suckers per arm. In total, 2,923,510 high-quality sequencing reads were obtained after filtering, with an average of 30,057 reads per paralarvae (ranging from 22 to 339,756 reads). Classification of the reads showed that 2,745,355 (93.9%) corresponded to O. vulgaris; 121,975 (4.17%) were identified as prey; and 56,180 (1.92%) were associated with contamination that was excluded from downstream analyses. Contaminants included fungi, humans, vertebrates, insects, as well as some marine species that were analysed in the laboratory in which the samples were processed, such as lobsters (see Appendix 1 for details). DNA of prey was detected in 95 out of 100 O. vulgaris paralarvae, with 1219 prey reads per paralarva on average (ranging from 1 to 29,677). The negative sample revealed 29 reads corresponding to O. vulgaris (Appendix 1).
The most abundant prey in terms of relative read abundance (RRA) across samples were the crabs Goneplax rhomboides (9.4%), Lophozozymus incisus (9.2%), Liocarcinus navigator (7.1%) and Liocarcinus corrugatus (6.8%) followed by the cnidarian Liriope tetraphylla (5.7%), the pteropod Cavolinia inflexa (5.5%) and the siphonophore Muggiaea atlantica (4.9%) (Appendix 2). The most frequent prey using FOO data were the crab G. rhomboides (detected in 26.3% of paralarvae), the siphonophore M. atlantica (24.2% of paralarvae), the crabs L. navigator and L. incisus (22.1% of paralarvae) and the cnidarian Clytia sp. (17.9% of paralarvae) (Fig. 1, Appendix 2 www.nature.com/scientificreports/ Diet along the Iberian Canary current upwelling system (ICC). Significant differences in the diet between the two sub-regions sampled within the ICC and also between the different locations sampled within www.nature.com/scientificreports/ those sub-regions, were identified using permutational multivariate analysis of variance (PERMANOVA; p < 0.001). Such diet differences can be observed in the PCO plots obtained with FOO ( Fig. 3a) and RRA ( Fig. 3b) data, where both factors ICC and location show similar patterns (covariation). The FOO plot shows great variability in the diets of the paralarvae and up to eleven species with correlations above 35% with PCO1 and PCO2 axes (Fig. 3a, blue vectors), with positive PCO1 values representing paralarvae collected in W Morocco and in the coastal location of NW Iberia, while negative values correspond with paralarvae collected from the ocean in NW Iberia. On the other hand, the RRA plot ( Fig. 3b) shows lower variability in the diets and eight species with correlations above 35%, but with similar distribution to that obtained with FOO data concerning the ICC subregions. The gradient from the coastal (positive values) to the oceanic areas (negative values) is more evident in this plot. Given that the diets are significantly different in the two ICC sub-regions, a detailed description is provided per region. ). An average of four prey (from one to 13) were detected in all octopus paralarvae except one sample with no prey DNA. Overall, 59 prey were identified: 35 crustaceans (16 crabs, seven copepods, three cladocerans, one hermit crab, three galatheids, one porcellanid, three euphausiids and one ostracod), seven fish, four siphonophores, four cnidarians, three molluscs (one pteropod, one cephalopod and one bivalve), two polychaetes, one ctenophore, one doliolid, one ophiuroid and a chaetognath (Fig. 1).
The diet of the paralarvae was significantly different according to the location collected (coast vs ocean, p = 0.001; factor correlated with PCO axis 1 in Fig. 3a,b) and the strata (p = 0.001). Pairwise analysis of the strata showed significant differences for the pairs < 5-< 100 m (p = 0.009) and < 5-< 500 m (p = 0.002). Pairwise analysis within the levels of the interaction factor (location and strata) revealed no significant differences.
The most discriminant prey that defined the two locations sampled off Northwest Iberian Peninsula (coast and ocean) are represented as pie charts in Fig. 4a. The siphonophore M. atlantica and the crabs L. navigator and G. rhomboides were the prey that better described the diets of coastal paralarvae. In contrast, the crabs L. incisus and L. corrugatus, and the copepod Centropages typicus were the better descriptors of the oceanic paralarvae (Fig. 4a). Eleven species are listed in the coastal area accounting for 90.7% of the variability observed and seven species in the ocean account for 91.7% of the total variability. The diets were different at the two locations  www.nature.com/scientificreports/ sampled, with only 18.3% and 18% of similarity in the coastal and oceanic areas, respectively (Fig. 4). The ten most discriminant organisms between the coastal and oceanic regions accounted for 50.7% of the total variability (Fig. 4b). The diet detected was markedly different as indicated by the 91.7% of dissimilarity observed between these two locations. DistLM results showed that 90% of the diet of the paralarvae collected along the Northwest Iberian Peninsula is predominantly represented by crabs (53.4%), followed by siphonophores (12.2%), copepods (12.3%), cnidarians (8.4%), and pteropods (3.7%) ( ). An average of 2.8 prey were detected from each paralarva (ranging from one to nine different prey), except for four samples with no prey reads. Overall, 45 prey were identified: 29 crustaceans (11 crabs, six euphausiids, seven copepods, one cladoceran and one ostracod), three fishes, five siphonophores, five cnidarians, two molluscs (one pteropod and one cephalopod) and one chaetognath (Fig. 1). The most frequent prey detected in the paralarvae were crabs (65.6% of prey reads, 18/31 paralarvae), cnidarians (2.7% of prey reads, 10/31 paralarvae), euphausiids and copepods (9.8 and 0.4% of prey reads, 9/31 paralarvae), and siphonophores, fish and non-octopus cephalopods (4.4, 0.2 and 0.2% of prey reads, 3/31 paralarvae) (Fig. 2b, Suppl. Table S1). The diet of the paralarvae collected in Morocco showed significant differences according to the location sampled (p = 0.049, Fig. 3b). Pairwise analysis of the factor location revealed significant differences in diets for the pair coast-ocean (p = 0.005). SIMPER analysis revealed that only three to five prey species account for more than 90% (between 90.1% to 97.1%) of the variability detected within the three different locations sampled in Morocco (Fig. 4a). The best descriptor species for the coast and coast-ocean locations were almost identical (the crabs L. marmoreus and G. rhomboides, the cnidarian L. tetraphylla and the euphausiid Nematoscelis megalops), with the oceanic location being more different (cyclopoid copepod, Brachioteuthis sp. and the siphonophore Rosacea sp.). The diets detected at the three locations gained diversity as the paralarvae drift offshore from the coastal area (24.4% similarity) within the upwelling filament and its periphery (C-O samples, 10% similarity) into the open ocean (4.7% similarity, Fig. 4a).
The importance of the top ten species that define the changes between locations is better observed when their relative contribution is ordered by their decreasing importance within the coastal group. The diets of the paralarvae at the three locations had very little in common (Fig. 4b), except for the crab L. marmoreus, the cnidarian L. tetraphylla and the euphausiid N. megalops, whose contribution decreased from the coastal into the oceanic locations. DistLM results showed that 95.6% of the diet of the paralarvae collected off Morocco is represented by crabs (46.2%), followed by copepods (23.1%), cnidarians (12.9%), euphausiids (9.3%), and fishes (4.2%) ( Table 1).
Ontogenetic changes in diet: trophic behaviour. Significant ontogenic differences between samples with a different number of suckers per arm (spa) was observed using FOO (%) data (PERMANOVA; p < 0.001). Pairwise analysis showed significant differences in the diet of paralarvae between 3 and 4 spa (p = 0.001), 3-> 5 (p = 0.004) and 4-> 5 suckers (p = 0.005). SIMPER analysis within the different groups of paralarvae revealed www.nature.com/scientificreports/ a decrease in prey diversity as the paralarvae grew. Pie charts in Fig. 5a show the species contributing > 90% of the total variability of the diet. Twelve species are responsible for 91.8% of the diet of paralarvae with three spa (collected from the coast and the ocean), whereas the diet of the paralarvae with four spa collected in the ocean is represented by eight prey (accounting for 92.5% of the variability), and six prey represent up to 94% of the diet of paralarvae with more than five spa. The ten most discriminant prey responsible for the differences detected between the diets of the three groups of suckers is shown in Table 2. The high dissimilarities obtained within each spa comparison (from 91.9 to 93.8%, Table 2) and the low similarities observed within each spa group (from 8.1 to 15.3%, Fig. 5a) suggest that the diets were very diverse. The importance of the top ten species that define the ontogenic changes is better observed when their relative contribution is ordered by their decreasing importance within the 3-spa group (Fig. 5b). It is evident that the contribution of meroplankton (i.e. larvae of benthic organisms like crabs) is replaced by holoplankton (copepods, siphonophores, cnidarians or pteropods) in older paralarvae. However, it is important to remember that the factors sucker number and location (coast, coast~ocean and ocean) are not independent. All paralarvae with 4 and > 5-spa were collected in the oceanic realm (that includes the coast-ocean and oceanic locations), while www.nature.com/scientificreports/ the paralarvae with 3-spa were the only ones collected in the coastal area but also in the two oceanic locations.
Since the diet in these locations was significantly different (p = 0.001; Figs. 3a, 4), the 3-spa group was split into coastal (3 spa C) and oceanic (3 spa O) paralarvae (Figs. 5c,d, 6). This way, it is possible to visualise the ontogenetic changes from the coastal area into the open ocean and identify the main prey groups driving those differences. Overall, meroplankton groups like bivalves, ophiuroids, hermit crabs, polychaetes, and porcellanid crabs were restricted to the coastal area in recently hatched paralarvae (3-spa). As they drift into the ocean, new holoplankton groups (copepods, siphonophores, cnidarians, euphausiids, pteropods, cladocerans and ostracods) progressively become more frequent in the diet (Fig. 5c,d).
The linear index of food selection (L) was calculated for every prey group ingested at the different areas sampled (Table 3), considering the FOO within the digestive tracts and the natural abundances of the prey groups ingested (Suppl . Table S2). This index provides an idea of prey group selectivity by comparing the relative contribution of the prey groups in the plankton and their presence in the digestive tracts. Positive values were obtained for 22 out of 26 prey groups evaluated (in green, Table 3), showing a preference for these groups. Negative values were obtained for ophiuroids, euphausiids, siphonophores and the copepod family Centropagidae, indicating ingestion of prey groups that were abundant in the zooplankton samples of Ría de Vigo (ophiuroids) and specific samples of NW Iberian Peninsula (euphausiids, Centropagidae and siphonophores) and W Morocco (Temoridae and Centropagidae).
The other selectivity index analyzed, called trophic niche breadth, allowed the evaluation of paralarval diets by comparing the frequency of distributions (based on FOO data) of prey in the diet and the relative proportion of the prey in the zooplankton. This measurement, derived using the Czekanowski's Index (CI), ranges from 0 to 1; in the context of diet, a value of 1 would indicate paralarvae target resources in proportion to their availability and reflects a generalist or non-specific predatory behaviour, whereas a value of 0 would indicate paralarvae are specialised exclusively on the rarest resource, and reflects a specialist predatory behaviour. The obtained CI ranged from 0.67 to 0.06, showing a broad range of trophic niches in the diets obtained for the paralarvae as they drift from the coast into the ocean (Fig. 6). Permanova analysis revealed no statistical differences in trophic niche breadth when paralarvae were grouped by sucker number (with averaged values of 0.25 ± 0.16 for 3-spa collected in the coastal area, 0.23 ± 0.12 for 3-spa collected in the ocean, 0.33 ± 0.17 for 4-spa and 0.24 ± 0.10 for > 5 suckers, Fig. 6 right corner). The averaged CI values suggest that O. vulgaris paralarvae can be considered specialist hunters with narrow trophic niches (CI < 0.5), feeding on prey groups that were not abundant in the plankton of the coastal and oceanic environments. However, within these groups, some paralarvae had broader trophic niches because they ingested prey that were abundant in specific samples (CI > 0.5), like euphausiids in sample S16 that accounted for 73% of the zooplankton abundance (Suppl. Table S2) or copepods of the family Centropagidae that dominated in most oceanic samples (Suppl . Table S2).

Discussion
In this study, we present the first insights into the diet of O. vulgaris paralarvae during their planktonic phase beyond the continental shelf in two different sub-regions of the ICC. Our data reveal marked differences in their diet within and between different sub-regions and ontogenic changes throughout development. Such changes in the diet are related to changes in feeding behaviour but also to changes in the zooplankton communities along coastal-ocean gradients. The sampling sites range from a coastal embayment including the Ría de Vigo 9 , to the continental shelf of the Northwest Iberian Peninsula 27-30 and the subtropical waters off Morocco [31][32][33][34] . These environments are markedly different in their physicochemical (temperature, salinity, transparency) and biological Table 1. DistLM results for the analysis of the main prey groups (FOO data) contributing to the diet of Octopus vulgaris paralarvae in the ICC zones sampled. The first two columns correspond to the marginal tests evaluating the independent contribution (Indep.) of the different prey groups. The following columns correspond to the sequential analysis of the groups retained by the model that maximised the % of explained variance (Adj. R 2 ). Prop. Indicates the individual % of variance explained by the groups as they are added into the model. %. Indicates the cumulative % of variance explained as prey groups are added into the model. Significant values are in bold. www.nature.com/scientificreports/ properties (primary production, zooplankton abundance and diversity), creating strong biological-and physical clines that paralarvae of O. vulgaris are subject to. The seasonal upwelling filament of Cape Silleiro (Northwest Iberian Peninsula, 42-43° N) and the permanent upwelling filament located off Cape Ghir, Morocco (30-31° N) export approximately 4 × 10 8 and 31 × 10 8 kg/year of organic carbon to the adjacent shelf representing 20 and 60% of the total phytoplankton primary production in these regions, respectively 35 . The upwelling filaments disrupt the strongly stratified oceanic waters creating distinct community gradients that contrast with the open ocean communities 36 . Octopus vulgaris paralarvae hatch in coastal areas and are transported by upwelling filaments hundreds of kilometres offshore to complete their planktonic phase 21 . As the paralarvae drift with the currents from the coast into the ocean, they are surrounded by zooplankton communities that gradually become more diverse but less abundant. In the oceanic environment beyond the continental shelf, 58 prey species were detected (39 being exclusive of this environment, Suppl. Table S1), including pteropods, ostracods, ctenophores, doliolids and oceanic cephalopods that were not detected before in Octopus gut contents.
Previous research suggests that the faunal composition in the ocean and upwelled waters of the ICC is dominated by small and medium-size copepods (accounting for 82.5-87% during the annual cycle) 36 . Despite  (Table 3), suggest that the paralarvae do not capture prey opportunistically based on their abundance (as expected for a generalist predator), but rather on other prey characteristics that can be related to their escape response or nutritional aspects. Copepods are fast swimmers with erratic movements that are difficult to predict and represent challenging prey for O. vulgaris. Copepod capture is a skill acquired throughout the planktonic stage, as observed with the squid Loligo opalescens in captivity 37 . Apart from prey characteristics, the anatomy of octopus paralarvae also limits their ability for prey hunting. Cephalopod paralarvae with a pair of feeding tentacles, such as loliginids or ommastrephids 1,2,11,38-40 , are better copepod hunters than O. vulgaris paralarvae which do not possess feeding tentacles. Herein, we observed that the paralarvae collected in the ocean consumed more copepods than those in the coastal area (Fig. 1, Suppl. Table 1). Since the paralarvae found in the open ocean are older (and have longer arms with an increased number of suckers) than those found near the coast 20 , this increase in copepod predation suggests that the paralarvae hunting skills increase as they develop, as it was observed for squid paralarvae 36 . Table 2. SIMPER analysis showing the FOO averaged abundance of the ten most discriminant species contributing to the ontogenic changes between paralarvae with different spa. The paralarvae were grouped into three groups (three suckers, 3 spa: n = 50; four suckers, 4 spa: n = 29; more than five suckers, > 5 spa: n = 16). Abbreviations: Diss/SD, discriminant power of each species; Var%, % of total variability explained by each species. www.nature.com/scientificreports/ The diets analysed showed that the diversity of prey is slightly higher in the oceanic area (57 different prey species detected in 62 paralarvae) than in the coastal area (51 different prey species detected in 33 paralarvae). However, this pattern is not consistent along the ICC: in the coastal and oceanic regions of the NW Iberian Peninsula, 46 and 28 prey species were detected in 25 and 39 paralarvae, respectively (Fig. 1). Conversely, 17 and 44 prey species were detected in eight and 23 paralarvae of Western Morocco's coastal and oceanic areas. The overall prey diversity detected in the coastal area is higher than previous studies, where 46 different prey taxa were detected in 56 out of 64 O. vulgaris with three suckers per arm with the COI gene 11 . Prey detected were very similar, except for the presence of polychaetes, cephalopods, doliolids and ctenophores, detected for the first time in this work. It is important to mention the absence of decapod families including Processidae, Alpheidae, Crangonidae and Thalassinidae, formerly detected in other dietary studies in O. vulgaris 8,11 but not present in this work. One possibility is that the modified primers, which included inosines to bind to any of the four nucleotides rather than degenerated bases 41,42 , may have prevented the amplification of those families.
The identification of discriminant species for the spatial and ontogenetic analyses carried out in this work was based on FOO rather than RRA. It is often assumed that because conversion to occurrence data moderates the impact of taxa-specific bias in marker signal, it provides a conservative view of diet. While it is true that occurrence-based summaries of diet are less affected by recovery bias and not skewed by read number, it has been observed that it is not necessarily true that they provide a more accurate representation of overall diet 43 . To evaluate this, we compared the discriminant species that defined the diets of the paralarvae at the different locations within the two areas samples of the ICC, using different datasets: FOO, RRA and log (x + 1) (Suppl. Fig. S1). FOO-based analysis recovered more discriminant species than log-based or RRA-based analyses, the latter being the most affected by read number. As an example, the coastal location of the NW Iberian Peninsula is defined by 10 species using FOO, seven species with log (x + 1) and five species with RRA. The primary drawback of occurrence data sets is that the importance of rare food taxa is often artificially inflated at the expense of food taxa eaten in large amounts 42 . However, there is a good concordance between the FOO and log (x + 1) analyses that identified similar discriminative power and explained variance for the top discriminant species. Contrarily, RRA-based analysis showed a reduced number of species with increased discriminatory power and explained variance due to the high read numbers of certain species.
The high diversity of euphausiids, cnidarians and siphonophores detected as prey in the ocean was remarkable, with siphonophores and cnidarians being the second and third most frequent prey ingested by O. vulgaris paralarvae (Fig. 2b). Natural contribution of gelatinous zooplankton like cnidarians and siphonophores is important inside the upwelling filaments, and their abundance gradually decreases as these mesoscale structures venture into the open ocean 44 . This trend was detected for the cnidarian Liriope tetraphylla in W Morocco, where a strong upwelling filament was sampled. Liriope was one of the top discriminant species defining the three locations sampled (Fig. 4a), with a relative contribution to the diet that decreased towards the oceanic samples (Fig. 4b).
In coastal areas, gelatinous zooplankton have also been detected as prey in O. vulgaris with three suckers 11 and loliginid paralarvae, like Alloteuthis media and Loligo vulgaris 11,38 . Different facts might help to explain both the predation and the importance of these organisms in cephalopod paralarvae. These organisms are proteinrich (18.3 ± 7.8% of dry mass 45 ), thus providing an easily digestible source of soluble nutrients 45,46 . Compared with highly motile prey like copepods, gelatinous plankton are relatively large slow drifters with predictable movements that are easy to capture by planktonic predators like cephalopod paralarvae 11,38 , fish larvae 47 and crustaceans 48,49 . According to the optimal foraging theory, the paralarvae will adopt a foraging strategy that provides the most benefit for the lowest cost, maximising the net energy gained. Within that scenario, a large gelatinous prey with lower energetic value would be easier to subdue than a more vigorous energy-rich prey (e.g. copepods), saving energy. Alternatively, considering the energy gained from the secondary copepod prey vs primary gelatinous prey, perhaps the octopus is specifically targeting the copepod rather than the gelatinous www.nature.com/scientificreports/ organism. Another important aspect of ingesting gelatinous prey is that transparency plays a key role in the open ocean for epipelagic organisms. This strategy is the most common solution to the dilemma of having nowhere to hide, a trait present in almost all marine phyla (e.g. salpids, crustaceans, leptocephalus, heteropods, polychaetes, chaetognaths or ctenophores, among others) inhabiting the epipelagic realm. Octopus paralarvae are completely transparent except for membranes enclosing the eyes and digestive gland, which are covered by reflective cells called iridophores that act as an ambient light reflector, concealing the opaque body organs 22 . The ingestion of transparent prey tend to make the paralarvae less conspicuous and more difficult to target by predators. This strategy of transparent predators ingesting transparent prey has also been suggested for the European eel larvae 50 and lobster larvae 45,48 . Crabs were the most common prey in all octopus analysed (68 out of 95 paralarvae analysed), even though their abundances in the ocean are an order of magnitude less (Table S2) than in shelf communities 9,10,33 . This crustacean group was preferentially ingested in all areas analysed (Table 3, green colour), and some of the species detected were among the most discriminant species to define both the spatial (Fig. 4a,b) and ontogenic variability in the diet (Fig. 5a,b). Interestingly, different crab larvae were consistently detected in paralarvae collected in www.nature.com/scientificreports/ both areas of the ICC, like Goneplax rhomboides or Liocarcinus navigator (Fig. 4a,b). These two species were mainly associated with coastal locations in both areas but also present within the coastal~ocean locations owed to the offshore transport within the upwelling filament. Most larvae cannot swim against horizontal currents but can swim faster than vertical currents to maintain a preferred depth, thus limiting their cross-shore dispersal 51 . Studies of larvae dispersal in the ICC suggest that vertical behaviour is one of the main biological mechanisms for crustacean larval retention over the shelf 28,29 , fish larvae retention 52 or fish larvae dispersal 53,54 . Upwelling filaments transport phytoplankton and zooplankton biomass between the productive coastal area and the oligotrophic oceanic realm 29,30,35 . Crab larvae are diel vertical migrators, and their horizontal distribution in the plankton (alongshore vs offshore) 28 greatly depends on their vertical distribution in the water 55  where Brachioteuthis riisei was one of the most abundant cephalopod paralarvae sampled 21 . The methodology applied in this work to characterise the diet prevents the detection of predation in conspecifics, but it may also constitute a potential prey since cannibalism has been detected in Octopus species 58,59 .
This study did not attempt to separate primary versus secondary prey items; items designated as prey could contribute either way. However, specific taxa detected in the digestive tracts of the paralarvae like rotifers, diatoms, or algae might result from secondary predation 60 . Considering Octopus paralarvae are visual predators that attack large prey compared with their size 22 , these groups are too small (rotifers or phytoplankton) or undigestible (algae) to be captured by the paralarvae. It is possible that these groups were ingested by other prey, like copepods or decapod larvae, shortly before these were captured by the paralarvae. The longer a prey item is inside a predator's gut, the harder it is to detect it because of DNA degradation. If the prey ingested by the paralarvae also ingested a prey item shortly before being captured, it could be possible to detect the DNA of both organisms inside the predator 61 . Molecular techniques are powerful tools to unravel trophic links, even in small organisms like O. vulgaris paralarvae. When dissections were carried out, only 9% of the paralarvae had amorphous contents inside the crop and stomach. However, prey DNA was successfully detected in 95%, highlighting the importance of molecular studies even when no prey contents are present.
Current aquaculture practices of octopus suffer from high levels of mortality, which may in part be due to sub-optimal diets in captivity 26 . These diets are based primarily on Artemia spp. (commonly used in aquaculture) enriched with commercial products or supplemented with other organisms like crustacean zoeae, copepods or amphipods. The increasing diversity of bacterial families detected in wild O. vulgaris paralarvae 24 , compared to the decreasing bacterial richness obtained in captive O. vulgaris, was hypothesised to be related to prey diversity available in the open ocean 24 . Our work shows that O. vulgaris paralarvae ingest numerous prey during the planktonic phase (at least 87 different taxa). Each prey has its own microbiome, thus adding complexity to the paralarvae's microbiome. This natural prey diversity in the wild heavily contrasts with the mono-diets commonly used in aquaculture, which offer less enriched nutrition and poor microbial communities, even though recently hatched octopus paralarvae possess a diverse microbial community 24 . However, diets based on Artemia profoundly impact the microbiota after a few days, leading to the expansion of opportunistic and pathogenic populations 24 . This has also been observed in other organisms reared in captivity, such as olive flounder 62 or white shrimp 63 . The O. vulgaris paralarvae studied in this work provide unique insight into their trophic ecology during the planktonic phase in the oceanic realm. This information can provide new avenues of research in captivity to reduce the high mortality levels currently constraining octopus aquaculture.

Material and methods
Oceanographic context. The Iberian-Canary current eastern boundary upwelling system (ICC) constitutes one of the world's oceans' four main eastern boundary upwelling systems. The ICC covers the latitudinal range between 12° and 43° N and can be broadly divided into five sub-regions according to their biogeographical characteristics 64 . This work is centred in three of these five sub-regions: the Galician, Portuguese and the Moroccan sub-regions (Fig. 7a). The NW Iberian Peninsula sub-region is the northernmost part of the ICC and includes the Galician and Portuguese sub-regions that are characterised by seasonal upwelling. During spring and summer (from March-April to September-October), north-easterly winds predominate in the Iberian basin, and mesoscale upwelling filaments develop intermittently (Fig. 7b) in association with irregularities in the coastline like capes [65][66][67] . On the other hand, the Moroccan sub-region experiences year-round upwelling that vary seasonally, with extended upwelling filaments (Fig. 7d), absence of freshwater inputs and massive dust inputs from the adjacent Sahara Desert 64,68 .
Zooplankton sampling. Plankton samples were collected during the multidisciplinary project "Canaries-Iberian Marine Ecosystem Exchanges (CAIBEX)" off the coast of Northwest Iberian Peninsula (CAIBEX-I, June 7th-24th) and Cape Ghir, West Morocco (CAIBEX-III, August 16th-September 5th) in 2009 onboard the research vessel (RV) "Sarmiento de Gamboa" (Fig. 7a). Two lagrangian experiments that followed two different water masses identified with drifting buoys were carried out off NW Iberian Peninsula: (1) an upwelling-relax- www.nature.com/scientificreports/ ation event in the open ocean in front of Portugal (L1: July 10th-13th inclusive Fig. 7b, samples S1-S7 in blue Fig. 7c); and (2) an upwelling event with alongshore transport along the continental shelf of Spain and Portugal, that can be identified by the sea surface temperature (SST) obtained at the end of the experiment with the violet colours representing the cold upwelled waters (L2: July 17th-20th inclusive, Fig. 7b, samples S13-S20 in green, Fig. 7c). Between L1 and L2, zooplankton samples were collected following a coastal-ocean gradient off the Portuguese coast (S8-S10) to observe the changes between these two environments, as well as two samples from the continental shelf of Galicia (S11 and S12). In contrast, a single Lagrangian experiment was carried out off the coast of W Morocco (31° N) following a strong upwelling filament 69,70 forced by the prevailing trade winds from the coast into the open ocean (L3: August 23rd-31st inclusive Fig. 7d, samples in orange, Fig. 7e). The magnitude of the cold upwelling filament was visible from the SST obtained from satellite data, as shown in Fig. 7d. Samples were also collected over the continental shelf (in green, Fig. 7e), in an area affected by the upwelled water over the continental slope (in orange, Fig. 7e), and in the open ocean (in blue, Fig. 7e). Mesozooplankton samples were collected day and night with two 750 mm diameter bongo nets equipped with 375 µm mesh and a mechanical flow meter. Three double-oblique towings were carried out at a ship speed www.nature.com/scientificreports/ of 2.5 knots over the continental slope (> 200 m depth): (1) from 0 to 500 m (< 500), (2) from 0 to 100 m (< 100), and (3) at the surface from 0 to 5 m (< 5). Over the continental shelf (< 200 m depth), only two double oblique towings were collected: from the surface down to 100 m (< 100 m, when the sea bottom was less than 100 m, the bongo net was lowered 10 m above the bottom) and at the surface (< 5). The bongo net was first lowered to the desired depth, towed for 30 min and subsequently hauled at 0.5 m s -1 . The net was recovered, cleaned on board, and placed back into the sea for the next towing. Plankton samples were fixed with 96% ethanol and stored at − 20 °C to facilitate DNA preservation. Mesozooplankton abundance was estimated for each sample by subsampling the original sample into an amount suitable for examination using a Folsom splitter. The subsample was made up to 300 ml, and several aliquots of 3 ml were obtained with a Stempel pipette, then identified and counted until at least 500 individuals were enumerated. Organisms were identified using a binocular (Nikon SMZ800) or inverted microscope (Nikon Eclipse TS100) to the lowest taxonomic level. All cephalopod paralarvae were separated from the zooplankton samples in the laboratory and stored individually in 70% ethanol at − 20 °C. In total, 134 O. vulgaris paralarvae were collected in NW Iberian Peninsula (n = 99) and W Morocco (n = 35), after which the number of suckers on the arms was counted. Of these, 60 paralarvae were chosen from the NW Iberian Peninsula (ranging from three to five suckers per arm) and 35 from Morocco (ranging from three to 15 suckers per arm) to study their diet (Table 4). Additionally, five recently hatched O. vulgaris paralarvae collected onboard the RV "Mytilus" at night at the surface of the Ría de Vigo (September 22nd 2010) were also included in the NW Iberian Peninsula group (n = 65). More detail of the sampling in the Ría de Vigo in 2010 can be found in 10 .
Library preparation and sequencing. The digestive tract of O. vulgaris paralarvae, which included the oesophagus, crop, stomach, caecum, digestive gland, and intestine, was dissected using entomological needles rinsed in ethanol and burned between every dissection. DNA was extracted using the QIAGEN DNeasy Blood and Tissue Kit, according to the manufacturer's instructions. A slight modification was made at the final elution stage of the extraction to increase the yield; the elution was repeated twice using two 20 µL aliquots of 45 °C ultrapure water and stored as a combined 40 µL eluate before use. A 300 bp region of the mitochondrial COI gene was amplified with NICO-F and NICO-R primers. The designed primers NICO-F and NICO-R (Suppl .  Table S3) are a modification of the degenerated primer mICOIintF 71 and HCO 72 , respectively. These primers were modified to add an inosine (I) that complements all four nucleotides 42 to capture a greater fraction of the prey. The primers contained a 5' tail containing a universal Illumina adapter sequence (italics) to enable multiplex indexing via primer extension PCR (Suppl . Table S3).
PCR reactions (12.5 μL) contained 10-30 ng of DNA template, 6.25 μL of RedTAQ ReadyMix (Sigma-Aldrich), 0.35 μL NICO-F (10 μM), 0.2 μL NICO-R (10 μM) and 0.1 μL MgCl 2 (25 nM). Touchdown-PCR was used to facilitate primer attachment with targets that were not 100% complementary. Cycling conditions consisted of an initial denaturation at 95 °C for 3 min followed by ten cycles of denaturation at 95 °C for 30 s, annealing at 57 °C for 30 s-decreasing 1 °C per cycle-and an extension at 72 °C for 40 s, then 25 cycles of 95 °C for 30 s, annealing at 47 °C for 30 s and an extension at 72 °C for 40 s, with a final step at 72 °C 4 min. PCR products were cleaned with AMPure beads according to the manufacturer's instructions and quantified with a Qubit HS kit in a Qubit™ fluorometer. Illumina multiplex indexes were added to the sequences in a second-round PCR before pooling 0.5 ng of each PCR product in a meta sample at a final concentration of 2.5 nM. To do so, some samples had to be diluted while others were repeated in duplicates and then concentrated to reach the same amount of DNA per sample. After adding 10% PhiX and diluting to 12.5 pM, the pooled library was sequenced using a v3 2 × 300 bp sequencing kit on an Illumina MiSeq. DNA extractions and PCR reactions were conducted in different rooms to avoid cross-contamination. Negatives were included for each round of PCR (with miliQ water instead of DNA) and in the meta-sample for sequencing. Quality filtering and bioinformatic analyses. Quality filtering was carried out following recommendations for Illumina platforms 73 . Reads that did not meet the following standards were removed: (1) Phred score below 30 (i.e., one error in 1000 bases), (2) less than 75% of target length, (3) less than three consecutive lowquality calls, (4) reads with ambiguous calls and (5) reads with less than five identical copies. The remaining paired-end reads were merged with PEAR v0.9.4 74 using a 99% sequence similarity threshold (minimum of five sequences per cluster). Merged reads were demultiplexed into individual sample read sets based on their corresponding adapter combination. Reads for which the indexes/primers did not match the expected sequences were www.nature.com/scientificreports/ discarded. The remaining reads were then filtered against a Kraken (v0.10.4) database to exclude archaeal and viral contamination 75 . The UCHIME algorithm of USEARCH (v 6.0.307) 76 was used to check for and remove chimeric sequences. Consensus sequences were BLASTed against sequences on GenBank, with the BLASTn function, and those with similarities above 80% were recovered. Amplicon sequencing variants (ASVs) corresponding to potential prey were assigned using the following criteria to taxonomical categories: ASVs with an identity higher than 97% were determined at the species level, ASVs between 93 and 97% were assigned to genus, ASVs between 93 and 90% were assigned to subfamily, and those below 90% were assigned to family or order 11 . Contamination reads were removed from downstream analyses to avoid confusion and included fungi, humans, vertebrates, insects, and some marine species analysed in the laboratory, like lobsters (see Appendix 1 for details).
Multivariate analysis of the diet. The software PRIMER6 and PERMANOVA+ 77 was used to analyse the diet of O. vulgaris paralarvae. Relative read abundance (RRA) and frequency of observance (FOO) were used as input data, and resemblance matrices were obtained with Bray-Curtis and Jaccard distances, respectively. An unbiased representation of the prey detected in the multidimensional space was visualised using principal coordinate analysis (PCO) plots. Permutational multivariate analysis of variance (PERMANOVA) tests using 9999 permutations were run to test the factors described in Table 4. Factor ICC with two levels: NW Iberian Peninsula, and W Morocco. Factor location with three levels: C, coast for all the samples collected over the continental shelf (green colour in Fig. 7c,e); C~O, for the samples collected over the continental slope but under the influence of the upwelling filament or coastal upwelled water (orange colour in Fig. 7e); and O, ocean for all the samples collected in the open ocean (blue colour in Fig. 7c,e). Factor strata with three levels; < 5 for samples collected between 0 and 5 m; < 100, between 0 and 100 m; and < 500, between 0 and 500 m. Factor day/night with two levels. Finally, the factor suckers, with three levels according to the number of spa: 3 (n = 51), 4 (n = 31), and > 5 (n = 18). PERMANOVA analysis was performed using the 'unrestricted permutation of raw data' option for these independent factors. Furthermore, a two-factor PERMANOVA test was also performed for the factors' location and strata. The factor strata was nested within the different locations, using the "permutation of residuals under the full model". The prey contributing most to similarities and dissimilarities between the paralarvae collected at different locations of the ICC and the paralarvae grouped by sucker count was determined using the program SIMPER 78 , using the Bray-Curtis dissimilarity matrix and FOO as input data to avoid biases related with read number. SIM-PER analysis allows the detection of the discriminant prey for the different factors analysed, their contribution to the total variability observed, and the discriminative power of the main prey driving the differences observed. Within each area of the ICC, the contribution of the different prey groups was analysed using FOO data with distance linear models (DistLM), applying an indicator that grouped the prey species at different taxonomic levels. A stepwise selection procedure was chosen using the adjusted R 2 as the selection criterion for retaining prey groups that summarised the diets of the different paralarvae at the different places of the ICC sampled.
Two trophic indices were analysed to evaluate the behavioural changes in trophic selection: the linear index of food selection (L) 79 for the different prey ingested and the trophic niche breadth (Czekanowski's index, CI) 80 for each paralarva. L ranges from − 1 to + 1, with positive values indicating a preference, negative values indicating avoidance or inaccessibility, and zero values showing random feeding 79 . This index was obtained for every prey detected with the formula (L = r i − p i ), where "r i " is the relative proportion of prey item i in the gut based on FOO and "p i " the relative proportion of prey item i in the zooplankton sample from which the paralarvae were selected. For this analysis, the mesozooplankton sample collected off the Ria de Vigo was considered separately from those of the NW Iberian Peninsula since they represent two zooplankton communities markedly different 9 . Prey were mainly grouped at the family level to unify the taxonomic levels obtained with the metagenomic study and the visual identification of the zooplankton. Furthermore, the trophic niche breadth for every octopus paralarvae was calculated using the CI 80 to evaluate the degree of overlap between the prey and their distribution in the zooplankton. Trophic niche breadth was calculated with the formula (CI = 1 − 0.5 ∑i | r i − p i |), where r i and p i are the relative abundances of resource item i eaten by the paralarvae based on FOO data (ri) and in the zooplankton (pi) (Suppl. Table S2). Values of CI range from 1 for the broadest possible niche (a population uses resources in proportion to their availability as a generalist predator) to [min p i ] for the narrowest possible niche (a population is specialised exclusively on the rarest resource). In the context of diet, we infer CI > 0.5 to reflect a generalist predator that targets prey in proportion to their availability, whereas CI < 0.5 reflects a specialist predator. Individual values of CI for each paralarvae were compared between the different ontogenetic groups (using spa as a factor) with PERMANOVA to test whether the foraging tactics of O. vulgaris changed with development.

Data availability
All data analysed during this study are included in this published article (and its Supplementary Information files). Raw files are available upon request.